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Abstract 

We compute the electric dipole moment of proton and neutron from lattice QCD simulations 
with Nf = 2 flavors of dynamical quarks at imaginary vacuum angle 9. The calculation proceeds 
via the CP odd form factor F-^. A novel feature of our calculation is that we use partially twisted 
boundary conditions to extract at zero momentum transfer. As a byproduct, we test the QCD 
vacuum at nonvanishing 9. 
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I. INTRODUCTION 



Current measurements of CP violating processes in the K and B meson sector would 
suggest that the phase of the CKM matrix provides a complete description. However, the 
baryon asymmetry of the universe cannot be described by this phase alone, suggesting that 
there are additional sources of CP violation awaiting discovery. 

QCD allows for a gauge invariant extra term in the action that is odd under CP trans- 
formations, 

S^S + i9Q, (1) 

where Q is the topological charge. Hence, there is the possibility of strong CP violation 
arising from a nonvanishing vacuum angle 9. The presence of CP violating forces implies a 
permanent electric dipole moment of the proton and neutron. This attribute is also deeply 
related to the question of baryon asymmetry of the universe 

The current experimental bound on the electric dipole moment of the neutron is 2| 

Ici^l < 2.9 X 10-^^efm (2) 

(Ifm = 10"^^ cm). Combining this bound with theoretical estimates of dj\f/9 allows us to 



le value of \9\. Current estimates from QCD sum rules jsl and 
J] give 1^1 ^ (1 — 3) X 10~^°. This anomaly is known as the 



derive an upper bound on t 
chiral perturbation theory 
strong CP problem. 

With the increasingly precise experimental efforts to observe the electric dipole moment 
of the neutron [sl, it is important to have a rigorous calculation directly from QCD. 

In this paper we present a calculation of in units of 9 with Nf = 2 flavors of dy- 
namical quarks using the lattice regularization. The novel feature of our work is that the 
simulations are jperformed directly at nonvanishing vacuum angle 9, in contrast to previous 

nnriQ □ 

lattice studies p, 17|, 18|, 19| (with the exception of [10]), which rely on reweighting correlation 



functions with topological charge that would otherwise vanish. The calculation becomes 



feasible if 9 is rotated to purely imaginary values [ll|, Il2|, Il3[. We expect from our method 
a much enhanced signal to noise ratio. In addition, we may hope to gain some insight into 
the dynamics of the 9 vacuum. For a recent review on this matter see jl^ . 
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II. THE ACTION 



The vacuum angle 9 can be rotated into the mass term, in the continuum and on the 



lattice 



15|, 



16|. This results in the fermionic action 

SF = i'{D + [cos{d /Nf) + I sm{9/Nf) 75] m} 7/^, (3) 

where summation over space-time coordinates {x, t} and quark flavors is understood, and 
D is the massless Dirac operator. For simplicity, we shall write in the following 

SF = i^{D + m + i{e/Nf)-f5m}il;, (4) 

with 

m = cos(6/Nf) m, 

(5) 

9 = tSin{9/Nf) Nf. 

We use clover fermions with Nf = 2 flavors of degenerate quarks. Taking into account 
that chiral symmetry is violated, we then have 

SF = ^{D + m + i {9r/2) Z^Zp-f5 m] ^, (6) 

where is the singlet renormalization constant of the vector Ward identity {VWI) quark 
mass and Zp that of the pseudoscalar density, and 9ji is the renormalized vacuum angle, 

9n={Z'^Zp)-' 9. (7) 

It can be shown that Zp is the same in both the singlet and nonsinglet case, while Zm is 
not Note that Z^Zp is scale independent, and in the continuum limit Z^Zp = 1 and 
9p^ = 9. 

III. THE SIMULATION 

The simulations are performed with the Iwasaki gauge action and said clover fermions 



with csw = 1-47 on 16^ 32 lattices at /5 = 2.1, k = 0.1357 18|]. The (massless) Dirac operator 



D in ([6]) is evaluated at Kc = 0.138984. The bare mass is taken to be am = 1/{2k) — 1/(2kc), 
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independent of 9. The resulting pion mass is m7r/mo ~ 0.8, corresponding to a quark mass 



oi m K, rris, mg being the strange quark mass. In 
was estimated at a ^ 0.11 fm, using the p mass to set the scale. 
The vacuum angle 6 is taken to be purely imaginary, 



181] the lattice spacing in the chiral limit 
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(8) 



resulting in the action 

SF = ij{D + m+{¥ /2)^^m]ij. (9) 

The simulations are done at 9^ = 0, 0.2, 0.4, 1.0 and 1.5, where we have collected 9000, 9000, 
7000, 6000 and 6000 trajectories of length one, respectively. 

We use the highly optimized HMC algorithm of the QCDSF Collaboration [3] for up- 
dating the gauge field. After integrating out the Grassmann fields, the action reads 



s[u, <p\ 0] = Sg[u] + Saot[u] + 0t(gtg)-V, 

where [t/] is the Iwasaki gauge action, (p^ and (j) are pseudofermion fields, and 



(10) 



5det [U] = -2 Tr log I 1 + Too + y 75 I , 



Q= ( i + r + y75 



Meo I l + r + y75 



Mo 



(11) 



with 



9' 



l--\9^ 

Kir 



(12) 



Meo and Moe are Wilson hopping matrices, which connect even with odd and odd with even 
sites, respectively, and T is the clover matrix 



T = -csw f^(^tivF^y{x). 



(13) 



We apply mass preconditioning a la Hasenbusch [20| an d sp lit the resulting action into 



three parts, each o 



which we put on separate time scales [2l|. We use Omelyan's second 



order integrator 22| to integrate Hamilton's equations of motion. 
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IV. RENORMALIZATION 



Let us now compute Z^Zp for our action and coupling, so that we can compare the 
results to phenomenology later on. We demand that the renormalized VWI and axial 
vector {AW I) quark masses are equal, 

niR = Z^m = ^ m = ffiR. (14) 

Zj 



p 



That gives Z^Zp = ZAfn/m. At the largest k value (smallest quark mass) of 18| we find 



m/m = 1.28(2). The renormalization constant Za has been computed nonperturbatively 



m 



231] for a variety of couplings. Extrapolating the numbers to 13 = 2.1 gives Z^ = 0.78(1). 
Multiplying these two pieces of information together, we obtain Z^Zp = 1.00(5). That 
means 6p = 6 to a. good precision. 

Alternatively, we may compute Z^Zp directly from Z^/Z^'^ and Z^^ Zp = Zp/Zg^ . As 
a comparison, the QCDSF Collaboration, using nonperturbatively improved clover fermions 



17| Z^/Z^'^ = 1.25(5) and |2| Zp/Z§ 



and the plaquette gauge action, finds at = 5.4 
0.81(2). Ahogether, this gives Z^Zp = 1.01(5), in agreement with the CP-PACS result, 
attesting clover fermions good chiral properties. 

V. CHARGE DISTRIBUTION AND 6 VACUUM 

Before we compute the electric dipole moment now, let us look at the distribution of 
topological charge and its dependence on 6. Having found that Oji = 6, the topological 
charge that follows from the chirally rotated action ([3]) is the so-called fermionic charge 

Q = mTr75M-\ (15) 



where M = D + m is the fermion matrix. The evaluation of ( 1151) requires the computation 



of the O(IOO) lowest 



demonstrated 



25 



26 



ying eigenvalues of M, which is numerically expensive. It has been 



271 ] that the fermionic charge and the so-called field theoretic charge, 
which is computed from the field strength tensor by applying an appropriate number of 
cooling sweeps and rounding the result to the nearest integer value, give consistent results. 
In the following we shall employ the field theoretic definition of the topological charge. It 
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Qh 0.06 




10 20 



FIG. 1: The topological charge distribution for 9^ = 0, 0.2, 0.4, 1.0 and 1.5, from right to left. To 
guide the eye, the distribution is compared to a Gaussian fit at = and 1.0. 



turns out that the results are rather independent of the degree of cooling. In fact, the 
numbers stabilize already after 0(10) cooling sweeps. The numbers quoted here refer to 
O(IOO) cooling sweeps. For a recent appraisal of the cooling method see 



In Fig. 1 we show the charge distribution for our five different values of 9^ . To identify the 
shape of the distribution, and to see how it changes with increasing value of r , we compare 
our data to a Gaussian fit at = and 1.0. A Gaussian distribution is the most common 



¥ 


(Q) 







-0.06(31) 


24.9(14) 


0.2 


-3.52(46) 


24.1(15) 


0.4 


-7.35(36) 


22.7(17) 


1.0 


-18.38(30) 


21.7(15) 


1.5 


-27.84(37) 


18.1(13) 



TABLE I: The average topological charge and charge squared for our values of 6^ . 
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distribution function for independent, randomly generated variables. 

In Table I we present the resulting average topological charge {Q) and charge squared, 

(Q2),^((g-(g)f) = (Q2)_(g)2^ (16) 

which we plot in Figs. 2 and 3. As with any distribution, besides the distributions mean, 
its skewness and kurtosis coefficients should be calculated in order to determine the type of 
distribution. We do so in Figs. 4 and 5, where we plot the skewness S, 

5 = |§^, {Q\^{{Q-{Q)r), (17) 

and kurtosis K, 

^^Wfc' (Q')^^{^Q-(Q))')-HiQ-(Q))"T- (18) 

Note that S and K have been normalized to {Q'^)c, different from the mathematical litera- 
ture, so that the volume dependence cancels out. 

Let us first look at the charge distributions in Fig. 1. At = the lattice data show a 
higher probability than a Gaussian distributed charge for intermediate values of IQI and a 
slightly thinner tail, while at larger values of 6^ the distributions appears to show a sharper 
peak and fatter tail. On top of that, the right tail becomes longer compared to the left one 
with increasing value of 6^ . 

This behavior is refiected in the skewness and kurtosis coefficients shown in Figs. 4 and 
5. Skewness is a measure of the degree of asymmetry of a distribution. If the right tail 
of the distribution is more pronounced than the left one, the distribution is said to have 
positive skewness, which is what we observe at > 0. If the reverse is true, it has negative 
skewness. Kurtosis is the degree of peakedness of a distribution. A distribution with positive 
kurtosis is called leptokurtic and has an acute peak around its mean. Examples of leptokurtic 
distributions include the Laplace distribution. A distribution with negative kurtosis is called 
platykurtic and has a smaller peak around its mean and a lower probability than a Gaussian 
distribution at extreme values. Such distributions are termed sub-Gaussian. ^ The kurtosis 
starts out negative at ^■'^ = and rises almost linearly to reach positive values at 6^ > 0.2. 



For recent quenched results see 



29 



30, 



3l| 



^ A distribution P{x) = exp[— /(x^)] is called sub-Gaussian if /'(x^) is strictly increasing on [0, oo). 
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FIG. 2: The average charge (O) compared to the reweighted numbers (x) and to the predictions 
of the Gaussian distribution (dotted hne) and the dilute instanton gas ()22p (dashed hne). 




FIG. 3: The average charge squared (O) compared to the reweighted numbers (x) and to the 
predictions of the Gaussian distribution (dotted hne) and the dilute instanton gas (j22|) (dashed 
line). 
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FIG. 4: The skewness compared to the predictions of the Gaussian distribution (dotted hne) and 
the dilute instanton gas (dashed hne). The errors shown are the naive ones, as our sample of large 
charges \Q — (Q) \ is too small to allow for a proper jackknive analysis, and may be underestimated. 



o 

A 
w 

V 



o 

A 
V 

II 




FIG. 5: The Kurtosis, together with its naive error, compared to the predictions of the Gaussian 
distribution (dotted line) and the dilute instanton gas (dashed line). 
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A popular model of the QCD vacuum is the dilute instanton gas [3^, whose charge 
distribution is given by a convolution of separate Poisson distributions for instantons and 
anti-instantons. Its free energy per spacetime volume V is 



F{9)=Xt{l-cos9), F{9) = -^\nZ{9), (19) 



where 



Xt 



That leads to 33 1 



V 



(20) 



(Q")c = -"^^, (21) 



which at imaginary 9 gives 



(22) 



{Q) = -Vxt smh9', 
{Q^)c = Vxt cosh^^ 
{Q^)c = -Vxt smh9', 
{Q^)c = Vxt cosh^^. 

In Figs. 2-5 we compare the lattice data for (Q), {Q'^)c, S and K with the predictions 
of the Gaussian distribution and the dilute instanton gas. While (Q) and {Q'^)c are in 
reasonable agreement with the results of the Gaussian and Poisson distribution for smaller 
values of 9^ , the higher cumulants S and K show a far different trend over the entire range 
of 9^ than the predictions of these simple models. 

Reweighting appears to be the accepted method for simulations at nonvanishing vacuum 
angle 9 and chemical potential fi. Having results of a direct simulation at nonvanishing value 
of 9 at hand, it is instructive to test how reliable the method actually is, given the fact that 
the simulations are necessarily restricted to a finite volume with limited absolute value of 
the topological charge. 

The reweighted charges are given by 

(^'^) = ^E^"^Q^"''^' E^« = l' (23) 

^(^) Q V 

where Pq denotes the probability of finding a configuration of charge Q in the ensemble of 
configurations, and 

Z{9) = J2Pq^~'"^- (24) 
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At imaginary 9 = —iO^ this becomes 

(^") = ^E^"^«^"''''' Z{e) = Y,PQe-'''^. (25) 

^ ' Q Q 

In Figs. 2 and 3 we compare (Q) and {Q'^)c with the reweighted numbers, where we have 
converted 6 to 6 using dS]). On the quantitative level, reweighting is not able to describe 
the data for 6^ > 0.4. The reason is that the reweighted charge distributions have largely 
the same form as the initial = distribution, but are merely shifted towards negative Q 
values, while the shape of the true charge distributions changes significantly with increasing 
value of 6^ . We may expect to find better agreement on larger volumes, provided Z{6) is 
analytic in 6. 

VI. NUCLEON FORM FACTORS AT 6* / 

At nonvanishing 6 the electromagnetic current between nucleon states can be decomposed 
in Euclidean space into 

{?', sV^ilP^ = U0{p', s') J^,ue{p, s), (26) 

with 



(27) 



where q = p' — p. The form factors and nucleon mass will generally depend on 6 with 
F^,^^ = F... and = ttln. The Dirac spinors are modified by a phase in the 6 vacuum. 



(28) 



ue{p,s)=e'''^'^^'uip,s), 
ue{p,s)=u{p,s)e'^^'^^\ 
so that the standard spinor relation is modified to 

As we are primarily interested in the electric dipole moment in the limit ^ — > 0, it is sufficient 
to consider the lowest order expansion only. Hence, we may write 

a{e) = a'e + o{e^). (30) 
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For our choice of 6 this then becomes to lowest order in 6 



^ue{p, s')ueip, s) 



2E 



N 



(31) 



Note that in Euchdean space = —{E' — E)'^ + [p' — p)'^ and 7p = ^-£'74 + 7p. 

We denote the two-point function of a nucleon of momentum p in the theta vacuum by 
G%j^{t,p). The phase factor a' of fl30l) can be obtained from the ratio of two-point functions 



Tr[G^^^(t;0)r4] 



-I Z,J2 p-rriNt 



2 ^—mivt 



(32) 



Tr[G'^^^(t;0)r475] - -c^'O' ^{Z^re 
where r4 = (1 -|- 74) /2. In Fig. 6 we show 

_ Tr[G^^^(t,0)r475] ^ 

for 6^ = 0.4. Fitting to a constant, we find 



(33) 







0.2 


0.048(3) 


0.4 


0.081(4) 



(34) 



The form factor F^lq"^), which is needed for the determination of the nucleon electric dipole 
moment, can be extracted from the ratio of three-point and two-point functions 

G%^jJt',t-p',p) 



Rf,{t',t;p',p) 



TT[G%^{t'-p')r4 



Tr [G%^{t;p')U Tr [G^^(t'; p')r4] Tr [G%^{t' - t;p)U 



Tr [G%j,it;p}T,] Tr [G%^it';p}T,] Tr [G%^it' - t;p')T,] 



1/2 



(35) 



EO'E^ 



where G^j^j^j^{{t',t]p',p) is the three-point function, with t' being the time location of the 
nucleon sink and t the time location of the current insertion, and the function F{r, J7^) is 



F(r,j-^) = -Trr 



E^''j4^ — Z7p' + m' 



X J. 



-E^74 — i'^p + m' 



N_ gia(e)75 



(36) 



12 



0.00 
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10 12 14 



FIG. 6: The ratio R{t) given in ([331) for 9^ = 0.4. 



with J'fj_ given in (p7j) . The three-point functions are calculated for various choices of nucleon 
polarization, F = F4, ^£47571 and ^F47572. The calculation follows the program of the 
QCDSF Collaboration 3J] for computing nucleon thr ee-p oint functions. We neglect finite 
volume corrections to two- and three-point functions 35|]. For we take the local vector 
current, which needs to be renormalized. We compute the corresponding renormalization 
from the proton form factor -Fi(O) at zero momentum transfer. 



constant Z 



V 



With conventional (periodic) boundary conditions momenta are quantized in units of 
2it/L, where L is the spatial extent of the lattice. For the lattices used in the current 
simulation, this means that the smallest nonvanishing momentum available is ~ 700 MeV. 
Since -F3 can only be computed at 7^ 0, we need to extrapolate to = to find d^. 
The momentum resolution of hadron observables can be significantly improved by varying 
the boundary conditions. It was demonstrated js^] that for processes without final state 
interactions, such as the form factors studied in this paper, it is sufficient to apply twisted 
boundary conditions to the valence quarks only. 

In our study we use partially twisted boundary conditions, i.e. combining gauge field 
configurations generated with sea quarks with periodic spatial boundary conditions with 
valence quarks with twisted boundary conditions. The boundary conditions of the valence 
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quarks attached to the electromagnetic current are 



k = 1,2,3. 



(37) 



By varying a we can tune the momenta of the nucleon continuously. We have chosen the 
following set of twist angles 

^ 2n 



a = —(0,0,0) 
a = — (0.36,0,0) 

Ij 

a = ^(0.36,0.36,0) 



(38) 



L 

2-K 



a 



(0.36,0.36,0.36) 



The dispersion relation for the nucleon then reads 

E = \i m 



N 



{p + dy 



We define the electric dipole moment as 



d 



eF|(0) 



N 



2m' 



(39) 



(40) 



N 



In Fig. 7 we show our results for for proton and neutron on configutaions with 6^ = 0.2 
(top) and = 0.4 (bottom). We find a very clean signal for the neutron form factor. The 
signal for the proton, on the other hand, is somewhat more noisy. The reason is that in this 
case one has to subtract -Ff(O) from the matrix element (12 7p . which is by far the largest 
contribution. 

To obtain the result at = 0, we first attempt a fit using a dipole ansatz 

^I(O) 



F (n 

which is indicated by the solid lines. At = we find 



(41) 



F3^(0)/(2m^) 



F«(0)/(2m;v) 



0.2 0.158(33) -0.108(17) (42) 
0.4 0.256(25) -0.193(12) 
Alternatively, if we assume that and Ff have similar behavior, then by forming the 
ratio 

(43) 



Ff'iq') 
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FIG. 8: The form factor ratio F^iq^) / Fi{q^) for proton and neutron, for 6^ = 0.2 (top) and 0^ = 0.4 
(bottom), respectively. 
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FIG. 9: The electric dipole moment for proton and neutron, together with a hnear plus cubic 
fit. The data points are horizontally displaced for better legibility. 



the renormalization constant Zy cancels, and we may hope to see a constant behavior as 
a function of g^. In Fig. 8 we show this ratio. Again, we find a very clean signal for the 
neutron. After performing a linear extrapolation to = 0, we obtain 



¥ 


Ff(0)/(2m^) 


F3«(0)/(2m;v) 


0.2 


0.141(16) 


-0.088(8) 


0.4 


0.257(18) 


-0.190(9) 



VII. ELECTRIC DIPOLE MOMENT 

Both values of ^3(0), (142!) and (l44l) . are consistent with each other within the error bars. 
In Fig. 9 we show the results together with a fit of the form 

'^'n = ^-^Q'^c¥\ (45) 
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as we are mainly interested in the derivative dd%/d6^ . The fit gives at = 

—Jf = 0.080(10) [e X fm] Proton, 

(46) 

= -0.049(5) [e X fm] Neutron. 

Ou 

Combining the upper experimental bound on the electric dipole moment of the neutron ([2]) 
with our result for dd%/d6^ , we may derive an upper bound on the vacuum angle 9. Taking 
our results at face value, we find 

1^1 < 6 X 10"^^ (47) 

It should be noted, however, that we are working at unphysically large quark mass yet, so 
that this result has limited phenomenological significance. 

VIII. CONCLUSION AND OUTLOOK 

We have performed simulations of QCD with A'^^ = 2 flavors of dynamical quarks at 
imaginary vacuum angle 6. It is the first time this has been done in full QCD. The use of 
partially twisted boundary conditions has allowed us to compute the proton and neutron 
form factor F^{q^) with high precision over the entire range of momenta down to (ag)^ ~ 0.02, 
which greatly facilitated the extrapolation to = 0. 

A further improvement of our calculation is that it does not require reweighting of the 
three-point functions (!35|) with the topological charge. Barring the fact that the lattice 
definition of topological charge is ambiguous, to some extent, reweighting does not describe 
the charge distribution accurately beyond 9^ > 0.2, which casts some doubts on the method, 
if taken at face value. 

Having demonstrated the benefit of simulations at imaginary 9, the next step is to extend 
the calculations to more realistic quark masses and larger lattices. The idea then is to make 
contact to the predictions of chiral perturbation theory, which allow for the extrapolation 



of the dipole moment from finite to infinite volume and to the physical quark mass |J, l38 | 



Chiral perturbation theory also predicts the 9 dependence of physical quantities, such as 



hadron masses. For the pion mass one finds 



39| 



ml{9) =ml{0) cos{9/Nf). (48) 
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FIG. 10: The pion mass squared as a function of 9 squared. The dashed hne shows the prediction 
of chiral perturbation theory. 



Though our quark mass is rather heavy, it is tempting to compare our results with fHHI) . 
This is done in Fig. 10. The variation of with 6 is found to be significant. Our results 
do not confirm the predictions of chiral perturbation theory. It will be interesting to see if 
this behavior persists at smaller quark masses. 

A caveat of our calculation is that clover fermions, though 0{a) improved, break chiral 
symmetry at finite lattice spacing. It is reassuring that Z:^Zp ^ 1, which indicates good 
chiral properties already. But there is the potential danger that the remaining O(a^) cor- 
rections will interfere wit^ 



rise to systematic errors 
smaller lattice spacing. ^ 



1 the assumed form of the nucleon matrix element fl35l) and give 
4ol |. To fully rule that out, we need to repeat the simulations at 



^ As an independent check for lattice artifacts, we have repeated the calculation for zero 6 angle of the 
valence quarks on our 9^ ~ 0.4 dynamical background field configurations and found a nonvanishing 
result for c?^. The result would have been zero in the absence of vacuum insertions of the pseudoscalar 



density 4fl |. 
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Finally, we plan to explore a far wider range of 6 values, intrigued by the results of a 
recent simulation of the 0(3) nonlinear sigma model in two dimension at imaginary 6 42|. 
By analytic continuation to real values of 6 it was possible to detect the phase transition of 
the model at 6' = vr and show that the mass gap vanishes at this point, in agreement with 
known results 43|. 
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